When not to avoid inbreeding: a gene’s eye view perspective
Author
Thomas Keaney, Arvid Agren and Hanna Kokko
Load packages
Code
# for tidy style coding and plottinglibrary(tidyverse) library(vroom) # to read lots of csv files at once# more table optionslibrary(pander) # for tableslibrary(kableExtra) # for scrolling tableslibrary(data.table) # for efficient handling of large dataframes# making ggplot more powerfullibrary(MetBrewer) # for colour palettes based upon artwork housed at the METlibrary(MoMAColors) # for colour palettes based upon artwork housed at MoMAlibrary(wesanderson) # for colour palettes based on wes anderson movieslibrary(rcartocolor) # for nice sequential colour schemeslibrary(PNWColors) # for colour palettes library(tidybayes) # for plotting distributionslibrary(stickylabeller) # labelling facets with strings in ggplotlibrary(geomtextpath) # for curved plot annotationslibrary(ggtext) # for markdown syntax in plot labelslibrary(patchwork) # for patching plots togetherlibrary(ggnewscale) # to reset scales in plots, allowing multiple fill arguments in ggplot# for computation speed checkslibrary(profvis) # breakdown of complex functionslibrary(bench) # individual functions
The model
To model the inclusive fitness consequences of inbreeding, four components that contribute to allele copy propagation need to be accounted for:
The number of offspring that inherit the focal allele when outcrossing: \(\frac{1}{2}(1 + F_\mathrm{A})n\), where \(F_\mathrm{A}\) is the inbreeding coefficient of the individual expressing the inbreeding gene.
The reduction in offspring produced due to inbreeding depression: \(\delta\)
The indirect fitness gain (offspring that inherit the focal allele from the focal individual’s mate) due to inbreeding: \(\frac{1}{2}(1 + F_\mathrm{R})(1-\delta)rn\), where \(r\) is the genealogical relatedness coefficient and \(F_\mathrm{R}\) is the inbreeding coefficient of the opposite-sex recipient individual that does not express the inbreeding gene.
the reproductive value \(v_i\) of each individual to the gene in question (\(i\) once again indexes the sex of each interactant).
A female’s inclusive fitness following a single inbred mating is:
, while fitness from a single outbred mating is simply \(v_\mathrm{f}\frac{1}{2}(1+F_\mathrm{f})n\). The ratio \(\frac{r(1+F_\mathrm{m})}{1+F_\mathrm{f}}\) is Hamilton’s (1972) gamete for gamete coefficient of relatedness, which when weighted by reproductive value \(\frac{v_\mathrm{m}r(1+F_\mathrm{m})}{v_\mathrm{m}(1+F_\mathrm{f})}\), is his life for life formulation of the relatedness coefficient (see Ubeda and Gardner 2017 for a similar use-case). In the autosomal case, \(v_\mathrm{f} = v_\mathrm{m}\) and there is no reason to a priori expect \(F_\mathrm{f} \neq F_\mathrm{m}\). This is not yet well expressed by my notation, but an additional, implicit assumption here is that the inbreeding coefficient for an inbreeding individual equals the inbreeding coefficient for an individual that outcrosses. For now let’s assume that is correct. I shall return to this point later. It follows that when \((1 + r)(1 - \delta) \gt 1\) the inclusive fitness of a female expressing the inbreeding allele exceeds mean female fitness, when considering only autosomal genes.
which for varying values of \(r\) looks like this:
Code
inbreeding_maximum_function <-function(r){ r / (1+ r)}parameters <-expand_grid(r =seq(from =0, to =1, by =0.05),delta =seq(from =0, to =1, by =0.05))r <- parameters %>%distinct(r)inbreeding_equilibria <-map_dfr(r, inbreeding_maximum_function) %>%rename(depression_threshold = r) %>%bind_cols(r)inbreeding_equilibria %>%ggplot(aes(x = r, y = depression_threshold)) +geom_line(linewidth =0.8) +coord_cartesian(ylim =c(0, 1)) +labs(x ='_r_, the individual-level relatedness coefficient',y =~delta~'(inbreeding depression)') +scale_x_continuous(expand =c(0, 0.009)) +scale_y_continuous(expand =c(0, 0)) +theme_bw() +theme(text =element_text(size =16),axis.title.x =element_markdown())
Code
# (prop fitness lost)\n that can be tolerated"
The parameter space above the curve shows where inbreeding avoidance is favoured by selection, while the parameter space below the curve shows where inbreeding preference is favoured.
Accounting for sex differences in genetic architecture
As stated above, inclusive fitness in the absence of inbreeding depression is \((1 + r)n\). Here, \(r\) represents the correlation between genotypes carried by interacting females and males, under the implicit assumption that both sexes have equal reproductive value for . That is, \(r\) is quantified at the level of the individual. However, given that there is sexual dimorphism in genetic architecture for many taxa, \(r\) does not sufficiently represent the correlation between genotypes for all loci. To model indirect allele propagation under inbreeding for specific regions of the genome, we multiply \(r\) with the new term \(a_{\mathrm{p},i}\), the probability that a locus present in the soma of one sex is also present in a gamete of an opposite sex mating partner of the \(i^{th}\) sex.
\(a_{\mathrm{p},i}\) quantifies the potential value of a mating partner for a focal allele. However, in the same way, we must also consider the value of the individual carrying the focal allele. Let this be \(a_{\mathrm{s},i}\), the probability that a locus present in the soma of of an individual is also present in a gamete of that individual. The direct component of fitness accrued becomes \(a_{\mathrm{s},i}n(1 - \delta)\), while indirect allele propagation is approximated by \(ra_{\mathrm{p},i}n(1 - \delta)\). Inbreeding is favoured in females when
where homozygosity is assumed to be high for outcrossing alleles in initially outbreeding populations.
To get a feel for what this modification does, consider an allele found at a diploid, autosomal locus. In this case, all of the gametes produced by the focal individual or a relative possess this locus, where they could potentially carry alleles identical by descent. Here, \(a_{\mathrm{s},i} = a_{\mathrm{p},i} = 1\) and the indirect component of inclusive fitness is dictated solely by \(r\). The results for this autosomal scenario are presented in Parker (1979), Kokko and Ots (2006) and others who have explored this topic. In contrast, an inbreeding preference allele present at a locus on a Y or W chromosome has no opportunity to propagate any alleles identical by descent through inbreeding, as these chromosomes are not carried by the gametes of the opposite sex mating partner. In this case \(a_{\mathrm{p},i} = 0\) when considering a relative’s gametes and \(a_{\mathrm{s},i} = 0.5\) with respect to the focal individual’s gametes. However, alleles on these hemizygous chromosomes still suffer the costs of inbreeding depression, which is a result of increased homozygosity for deleterious recessive alleles throughout the genome. Conflict over the expression of inbreeding preference between alleles present on autosomes and those present on hemizygous sex chromosomes is immediately clear.
X- or Z-linked loci present an interesting intermediate case, with sex-specific values for \(a_{\mathrm{p},i}\) and \(a_{\mathrm{s},i}\). When the inbreeding locus is carried by the sex with homozygous sex chromosomes, \(a_{\mathrm{s},i} = 1\), but \(a_{\mathrm{p},i}\) is half that of autosomal loci, whereas when the locus is found within the hemizygous sex \(a_{\mathrm{s},i} = 0.5\) and \(a_{\mathrm{p},i} = 1\). Using loci on the X as an example, those present in a XX female are only found in ~50% of a interacting males gametes, as the remaining 50% carry Y chromosomes (assuming an even primary sex ratio). When an X-linked locus is found in a male, an interacting female’s gametes all carry X chromosomes and \(a_{\mathrm{p,f}} = 1\).
The X/Z situation appears additionally complex because there may be an element of frequency dependence to the kin selected benefits. When an inbreeding allele on an autosome is rare, then the chance of a relative carrying two copies is low, whereas when the allele is common, this chance is much higher. Rarity could therefore lead to similar fitness outcomes for autosomal and X/Z linked alleles (when present in the hemizygous sex), while commonality of the allele likely roughly equates to the conflicting situation outlined in the above paragraph. However, alleles favoured by kin selection rapidly become common within families, even whilst rare across the population (Fisher, 1930). This effect is especially pronounced with inbreeding, which independently increases the frequency of homozygotes within families. In Section 3, we simulated this for autosomal, X- and Z-linked loci which confirmed that homozygotes rapidly become the most common genotype within inbreeding families.
Table 1. Values of the parameter \(a_{\mathrm{s},i}\) and \(a_{\mathrm{p},i}\) parameters for different regions of the genome. Note that cytoplasmic chromosomes are assumed to have exclusive maternal inheritance.
Code
tibble(`Prob. that focal individual gametes carry focal locus`=c(1, 1, 1, 0.5, 0.5, 0),`Prob. that opposite sex gametes carry focal locus`=c(1, 0.5, 0, 1, 0, 1),`Relevant cases`=c("Autosomes in either sex", "X chromosomes in females, Z chromosomes in males","Cytoplasmic chromosomes in females","X chromosomes in males, Z chromosomes in females","Y chromosomes in males, W chromosomes in females","Cytoplasmic chromosomes in males")) %>%pander(split.cell =20, split.table =Inf)
Prob. that focal individual gametes carry focal locus
Prob. that opposite sex gametes carry focal locus
Relevant cases
1
1
Autosomes in either sex
1
0.5
X chromosomes in females, Z chromosomes in males
1
0
Cytoplasmic chromosomes in females
0.5
1
X chromosomes in males, Z chromosomes in females
0.5
0
Y chromosomes in males, W chromosomes in females
0
1
Cytoplasmic chromosomes in males
Once again we can find the condition where breeding with a relative returns greater fitness than an inbreeding avoidance strategy, this time accounting for genetic architecture
Figure XX. curves show the extent of inbreeding depression, expressed as a reduction in offspring produced, at which inbreeding no longer increases inclusive fitness. \(~\)
Differences between the sexes beyond genetic architecture
In his 1979 book chapter, Parker considered the inclusive fitness results of breeding with a relative and identified that females and males should have different tolerances for inbreeding depression. The key departure from the unlimited polygyny case presented in the above equations is that a cost to future reproductive success is included for males, i.e. due to finite sperm production, parental care, or harmful mating behaviour such as sexual cannibalism.
Parker modelled the condition for monandrous females to prefer incestuous matings when also presented with an outcrossing opportunityas
\[n(1 - \delta) + rn(1-\delta) - crn \gt n\] where the new third term on the lefthand side of the equation is the number of alleles that were not indirectly propagated by her male relative by future outcrossing. \(c\) is the cost of the present mating, relative to what is lost for a female. This can be considered a ratio of parental investment. When \(c = 1\) parental investment in the current bout of reproduction is even between the sexes. Alternatively, if males only contribute cheaply produced sperm to an incestuous mating, the cost of mating is likely very small relative to females i.e. \(c ~ 0\).
We add \(a_\mathrm{f}\) and \(a_\mathrm{m}\) to the equation and letting \(n = 1\), reduce to
\[\delta_\mathrm{m} = \frac{a_\mathrm{s,m}}{a_\mathrm{s,m} + ra_\mathrm{p,f}},\] the inbreeding depression threshold when mating is cheap for males.
Plotting sexual and intragenomic conflicts
When is inbreeding favoured in each sex?
In the figure below, the curves indicate the level of inbreeding depression that can be tolerated for a given value of individual-level relatedness. The plot is split into panels by \(c\), the cost of mating for males relative to females.
To quantify conflict in its various forms we use a version of the I index presented in Innocenti and Morrow (2010):
\[I = \frac{s_\mathrm{i} s_\mathrm{j}}{\sqrt{\frac{{s_\mathrm{i}}^2 + {s_\mathrm{j}}^2}{2}}}\] where \(s_\mathrm{i}\) and \(s_\mathrm{j}\) are the selection coefficients that we calculate for inbreeding alleles in the \(i^{th}\) and \(j^{th}\) context, respectively. \(I\) is most extreme when selection is strong in both contexts, positive when selection is context concordant and negative when selection is antagonistic.
Let’s plot the conflict zones between chromosomes, within the sexes:
Code
oranges <-c("#ffe6b7", "#ffd06f", "#f7aa58", "#ef8a47", "#e76254") Intragenomic_conflict_females <- analytical_results %>%select(1:5, female_fitness_contrast) %>%mutate(Case =case_when(a_m ==0& a_f ==1~"Cytoplasmic", a_m ==0& a_f ==0.5~"W", a_m ==0.5& a_f ==1~"X", a_m ==1& a_f ==0.5~"Z", a_m ==1& a_f ==1~"Autosome")) %>%filter(!is.na(Case)) %>%select(-c(a_f, a_m)) %>%pivot_wider(names_from = Case, values_from = female_fitness_contrast) %>%mutate(`Autosome vs X`=I_calculator(Autosome, X),`Z vs Autosome`=I_calculator(Autosome, Z),`Autosome vs W`=I_calculator(Autosome, `W`),`Autosome vs Cytoplasmic`=I_calculator(Autosome, `Cytoplasmic`),`Z vs W`=I_calculator(Z, `W`),`Z vs Cytoplasmic`=I_calculator(Z, `Cytoplasmic`),`X vs Cytoplasmic`=I_calculator(X, `Cytoplasmic`),`W vs Cytoplasmic`=I_calculator(W, `Cytoplasmic`)) %>%pivot_longer(cols =9:16, names_to ="contrast", values_to ="Evolutionary_conflict") %>%mutate(relationship =case_when( contrast =="Autosome vs X"& Autosome >0& X >0~"Inbreeding favoured in both contexts", contrast =="Autosome vs X"& Autosome <0& X <0~"Inbreeding deleterious in both contexts", contrast =="Z vs Autosome"& Autosome >0& Z >0~"Inbreeding favoured in both contexts", contrast =="Z vs Autosome"& Autosome <0& Z <0~"Inbreeding deleterious in both contexts", contrast =="Autosome vs W"& Autosome >0&`W`>0~"Inbreeding favoured in both contexts", contrast =="Autosome vs W"& Autosome <0&`W`<0~"Inbreeding deleterious in both contexts", contrast =="Autosome vs Cytoplasmic"& Autosome >0&`Cytoplasmic`>0~"Inbreeding favoured in both contexts", contrast =="Autosome vs Cytoplasmic"& Autosome <0&`Cytoplasmic`<0~"Inbreeding deleterious in both contexts", contrast =="Z vs W"& Z >0&`W`>0~"Inbreeding favoured in both contexts", contrast =="Z vs W"& Z <0&`W`<0~"Inbreeding deleterious in both contexts", contrast =="Z vs Cytoplasmic"& Z >0&`Cytoplasmic`>0~"Inbreeding favoured in both contexts", contrast =="Z vs Cytoplasmic"& Z <0&`Cytoplasmic`<0~"Inbreeding deleterious in both contexts", contrast =="X vs Cytoplasmic"& X >0&`Cytoplasmic`>0~"Inbreeding favoured in both contexts", contrast =="X vs Cytoplasmic"& X <0&`Cytoplasmic`<0~"Inbreeding deleterious in both contexts", contrast =="W vs Cytoplasmic"& W >0&`Cytoplasmic`>0~"Inbreeding favoured in both contexts", contrast =="W vs Cytoplasmic"& W <0&`Cytoplasmic`<0~"Inbreeding deleterious in both contexts",.default ="Intragenomic conflict")) %>%mutate(contrast =fct_relevel(contrast, "W vs Cytoplasmic","Autosome vs X", "Z vs Autosome","X vs Cytoplasmic","Autosome vs W", "Autosome vs Cytoplasmic","Z vs W","Z vs Cytoplasmic")) %>%filter(c ==0) # remove if we want more c valuesIntragenomic_conflict_males <- analytical_results %>%select(1:5, male_fitness_contrast) %>%mutate(Case =case_when(a_m ==0& a_f ==1~"Cytoplasmic", a_m ==0.5& a_f ==1~"X", a_m ==0.5& a_f ==0~"Y", a_m ==1& a_f ==0.5~"Z", a_m ==1& a_f ==1~"Autosome")) %>%filter(!is.na(Case)) %>%select(-c(a_f, a_m)) %>%pivot_wider(names_from = Case, values_from = male_fitness_contrast) %>%mutate(`Autosome vs X`=I_calculator(Autosome, X),`Z vs Autosome`=I_calculator(Autosome, Z),`Y vs Autosome`=I_calculator(Autosome, Y),`Autosome vs Cytoplasmic`=I_calculator(Autosome, Cytoplasmic),`Y vs X`=I_calculator(Y, X),`X vs Cytoplasmic`=I_calculator(X, Cytoplasmic),`Z vs Cytoplasmic`=I_calculator(Z, Cytoplasmic),`Y vs Cytoplasmic`=I_calculator(Y, Cytoplasmic),) %>%pivot_longer(cols =9:16, names_to ="contrast", values_to ="Evolutionary_conflict") %>%mutate(relationship =case_when( contrast =="Autosome vs X"& Autosome >0& X >0~"Inbreeding favoured in both contexts", contrast =="Autosome vs X"& Autosome <0& X <0~"Inbreeding deleterious in both contexts", contrast =="Z vs Autosome"& Autosome >0& Z >0~"Inbreeding favoured in both contexts", contrast =="Z vs Autosome"& Autosome <0& Z <0~"Inbreeding deleterious in both contexts", contrast =="Y vs Autosome"& Autosome >0& Y >0~"Inbreeding favoured in both contexts", contrast =="Y vs Autosome"& Autosome <0& Y <0~"Inbreeding deleterious in both contexts", contrast =="Autosome vs Cytoplasmic"& Autosome >0&`Cytoplasmic`>0~"Inbreeding favoured in both contexts", contrast =="Autosome vs Cytoplasmic"& Autosome <0&`Cytoplasmic`<0~"Inbreeding deleterious in both contexts", contrast =="Y vs X"& X >0& Y >0~"Inbreeding favoured in both contexts", contrast =="Y vs X"& X <0& Y <0~"Inbreeding deleterious in both contexts", contrast =="X vs Cytoplasmic"& X >0&`Cytoplasmic`>0~"Inbreeding favoured in both contexts", contrast =="X vs Cytoplasmic"& X <0&`Cytoplasmic`<0~"Inbreeding deleterious in both contexts", contrast =="Z vs Cytoplasmic"& Z >0&`Cytoplasmic`>0~"Inbreeding favoured in both contexts", contrast =="Z vs Cytoplasmic"& Z <0&`Cytoplasmic`<0~"Inbreeding deleterious in both contexts", contrast =="Y vs Cytoplasmic"& Y >0&`Cytoplasmic`>0~"Inbreeding favoured in both contexts", contrast =="Y vs Cytoplasmic"& Y <0&`Cytoplasmic`<0~"Inbreeding deleterious in both contexts",.default ="Intragenomic conflict")) %>%mutate(contrast =fct_relevel(contrast, "Autosome vs X", "Z vs Autosome","Y vs Autosome","Y vs X", "X vs Cytoplasmic","Autosome vs Cytoplasmic","Z vs Cytoplasmic", "Y vs Cytoplasmic")) %>%filter(c ==0) # remove if we want more c valuesmake_genomic_conflict_plot <-function(data, enter_title, colour_pal){ data %>%ggplot(aes(x = r, y = D)) +geom_blank() +geom_tile(data = data %>%filter(relationship =="Intragenomic conflict"),aes(fill = Evolutionary_conflict*-1)) +#geom_tile(data = data,# aes(fill = Evolutionary_conflict)) + scale_fill_gradientn(colours = colour_pal, limits =c(0, 0.6), #na.value = "white",breaks =c(0, 0.6),labels =c("Weaker conflict", "Stronger conflict")) +labs(fill ="Evolutionary concordance") +new_scale_fill() +geom_tile(data = data %>%filter(relationship !="Intragenomic conflict"),aes(fill = relationship), alpha =0.75) +scale_fill_manual(values =c("#fbe6c5", "#d2fbd4"), labels =c("Inbreeding deleterious\nin both contexts", "Inbreeding favoured\nin both contexts")) +stat_contour(aes(z = Evolutionary_conflict*-1), colour ="black", binwidth =25,breaks =c(0, 0.2, 0.4, 0.6)) +facet_wrap(~contrast, nrow =3,scales ="free", strip.position =c("top"),labeller =label_glue('{`contrast`}')) +labs(fill ="Evolutionary concordance",y =~delta~'(inbreeding depression)', x ='_r_, the individual-level relatedness coefficient',title = enter_title) +scale_x_continuous(expand =c(0, 0), breaks =c(0, 0.25, 0.5, 0.75)) +scale_y_continuous(expand =c(0, 0)) +# labels = c(0, 25, 50, 75, 90)) +theme(panel.border =element_rect(fill =NA, colour ="black", size = .8),panel.grid.minor =element_blank(),strip.background =element_rect(colour ="black", fill ="Aliceblue", linewidth = .8),axis.title.x =element_markdown(),plot.title =element_text(hjust =0.5, size =12)) }icf <-make_genomic_conflict_plot(Intragenomic_conflict_females, "A. Female inter-chromosomal conflict", oranges)icm <-make_genomic_conflict_plot(Intragenomic_conflict_males, "B. Male inter-chromosomal conflict", oranges)
Code
icf
Figure XX. the parameter space where the expression of inbreeding in females increases the propagation of one type of chromosome, but has an opposite, deleterious affect on the propagation of a second chromosome class; a so called inter-chromosomal conflict.
The plot shows regions of inter-chromosomal conflict within females, where chromosome classes with a high probability of encountering opposite-sex gametes carrying the inbreeding locus have high fitness in the conflict zone. In some cases, the fitness interests of chromosome classes align e.g. cytoplasmic and W chromosomes.
Code
icm
Figure XX. the parameter space where the expression of inbreeding in males increases the propagation of one type of chromosome, but has an opposite, deleterious affect on the propagation of a second chromosome class.
Here inter-chromosomal conflict is plotted when males express the inbreeding allele. Importantly, the direction of selection for inbreeding is reversed in comparison to the female case in the zones of conflict. In males, chromosome classes with high probabilities of inbreeding loci being found in female gametes now have low fitness in the conflict space.
Now let’s have a look at chromosome specific conflict between the sexes
Figure XX. conditions where loci are expected to be under intra-chromosomal sexual conflict over inbreeding preference. In regions of sexual conflict, inbreeding preference is always favoured in males, but has negative fitness consequences if expressed by females (assuming that males invest less into mating than females). Note that intra-chromosomal conflict encompasses both the intra- and inter-locus forms of sexual conflict.
Make Figure SX - fix this then analytics done
Also equation numbering
Code
# calculate every combinationIntragenomic_all_combos <- analytical_results %>%select(1:5, female_fitness_contrast, male_fitness_contrast) %>%mutate(Case =case_when(a_m ==0& a_f ==0.5~"W", a_m ==0& a_f ==1~"Cytoplasmic", a_m ==0.5& a_f ==1~"X", a_m ==0.5& a_f ==0~"Y", a_m ==1& a_f ==0.5~"Z", a_m ==1& a_f ==1~"Autosome")) %>%filter(!is.na(Case) & c ==0) %>%select(-c(a_f, a_m)) %>%pivot_wider(names_from = Case, values_from =c(female_fitness_contrast, male_fitness_contrast)) %>%mutate(`Autosome (female) vs Autosome (male)`=I_calculator(female_fitness_contrast_Autosome, male_fitness_contrast_Autosome),`Autosome (female) vs X (female)`=I_calculator(female_fitness_contrast_Autosome, female_fitness_contrast_X),`Autosome (female) vs X (male)`=I_calculator(female_fitness_contrast_Autosome, male_fitness_contrast_X),`Autosome (female) vs Y (male)`=I_calculator(female_fitness_contrast_Autosome, male_fitness_contrast_Y),`Autosome (female) vs Z (female)`=I_calculator(female_fitness_contrast_Autosome, female_fitness_contrast_Z),`Autosome (female) vs Z (male)`=I_calculator(female_fitness_contrast_Autosome, male_fitness_contrast_Z),`Autosome (female) vs W or Cytoplasmic (female)`=I_calculator(female_fitness_contrast_Autosome,`female_fitness_contrast_W or Cytoplasmic`),`Autosome (female) vs Cytoplasmic (male)`=I_calculator(female_fitness_contrast_Autosome,`male_fitness_contrast_W or Cytoplasmic`),`X (female) vs Autosome (male)`=I_calculator(female_fitness_contrast_X, male_fitness_contrast_Autosome),`X (female) vs X (male)`=I_calculator(female_fitness_contrast_X, male_fitness_contrast_X),`X (female) vs Y (male)`=I_calculator(female_fitness_contrast_X, male_fitness_contrast_Y),`X (female) vs Z (female)`=I_calculator(female_fitness_contrast_X, female_fitness_contrast_Z),`X (female) vs Z (male)`=I_calculator(female_fitness_contrast_X, male_fitness_contrast_Z),`X (female) vs W or Cytoplasmic (female)`=I_calculator(female_fitness_contrast_X,`female_fitness_contrast_W or Cytoplasmic`),`X (female) vs Cytoplasmic (male)`=I_calculator(female_fitness_contrast_X,`male_fitness_contrast_W or Cytoplasmic`),`Z (female) vs Autosome (male)`=I_calculator(female_fitness_contrast_Z, male_fitness_contrast_Autosome),`Z (female) vs X (male)`=I_calculator(female_fitness_contrast_Z, male_fitness_contrast_X),`Z (female) vs Y (male)`=I_calculator(female_fitness_contrast_Z, male_fitness_contrast_Y),`Z (female) vs Z (male)`=I_calculator(female_fitness_contrast_Z, male_fitness_contrast_Z),`Z (female) vs W or Cytoplasmic (female)`=I_calculator(female_fitness_contrast_Z, `female_fitness_contrast_W or Cytoplasmic`),`Z (female) vs Cytoplasmic (male)`=I_calculator(female_fitness_contrast_Z, `male_fitness_contrast_W or Cytoplasmic`),`W or Cytoplasmic (female) vs Autosome (male)`=I_calculator(`female_fitness_contrast_W or Cytoplasmic`, male_fitness_contrast_Autosome),`W or Cytoplasmic (female) vs X (male)`=I_calculator(`female_fitness_contrast_W or Cytoplasmic`, male_fitness_contrast_X),`W or Cytoplasmic (female) vs Y (male)`=I_calculator(`female_fitness_contrast_W or Cytoplasmic`, male_fitness_contrast_Y),`W or Cytoplasmic (female) vs Z (male)`=I_calculator(`female_fitness_contrast_W or Cytoplasmic`, male_fitness_contrast_Z),`W or Cytoplasmic (female) vs Cytoplasmic (male)`=I_calculator(`female_fitness_contrast_W or Cytoplasmic`,`male_fitness_contrast_W or Cytoplasmic`),`Autosome (male) vs X (male)`=I_calculator(male_fitness_contrast_Autosome, male_fitness_contrast_X),`Autosome (male) vs Y (male)`=I_calculator(male_fitness_contrast_Autosome, male_fitness_contrast_Y),`Autosome (male) vs Z (male)`=I_calculator(male_fitness_contrast_Autosome, male_fitness_contrast_Z),`Autosome (male) vs Cytoplasmic (male)`=I_calculator(male_fitness_contrast_Autosome, `male_fitness_contrast_W or Cytoplasmic`),`X (male) vs Y (male)`=I_calculator(male_fitness_contrast_X, male_fitness_contrast_Y),`X (male) vs Z (male)`=I_calculator(male_fitness_contrast_X, male_fitness_contrast_Z),`X (male) vs Cytoplasmic (male)`=I_calculator(male_fitness_contrast_X, `male_fitness_contrast_W or Cytoplasmic`),`Y (male) vs Z (male)`=I_calculator(male_fitness_contrast_Y, male_fitness_contrast_Z),`Y (male) vs Cytoplasmic (male)`=I_calculator(male_fitness_contrast_Y, `male_fitness_contrast_W or Cytoplasmic`),`Z (male) vs Cytoplasmic (male)`=I_calculator(male_fitness_contrast_Z, `male_fitness_contrast_W or Cytoplasmic`) ) %>%pivot_longer(cols =14:49, names_to ="contrast", values_to ="Evolutionary_conflict") %>%mutate(relationship =case_when(# beneficial in both contexts contrast =="Autosome (female) vs Autosome (male)"& female_fitness_contrast_Autosome >0& male_fitness_contrast_Autosome >0| contrast =="Autosome (female) vs X (female)"& female_fitness_contrast_Autosome >0& female_fitness_contrast_X >0| contrast =="Autosome (female) vs X (male)"& female_fitness_contrast_Autosome >0& male_fitness_contrast_X >0| contrast =="Autosome (female) vs Y (male)"& female_fitness_contrast_Autosome >0& male_fitness_contrast_Y >0| contrast =="Autosome (female) vs Z (female)"& female_fitness_contrast_Autosome >0& female_fitness_contrast_Z >0| contrast =="Autosome (female) vs Z (male)"& female_fitness_contrast_Autosome >0& male_fitness_contrast_Z >0| contrast =="Autosome (female) vs W or Cytoplasmic (female)"& female_fitness_contrast_Autosome >0&`female_fitness_contrast_W or Cytoplasmic`>0| contrast =="Autosome (female) vs Cytoplasmic (male)"& female_fitness_contrast_Autosome >0&`male_fitness_contrast_W or Cytoplasmic`>0| contrast =="X (female) vs Autosome (male)"& female_fitness_contrast_X >0& male_fitness_contrast_Autosome >0| contrast =="X (female) vs X (male)"& female_fitness_contrast_X >0& male_fitness_contrast_X >0| contrast =="X (female) vs Y (male)"& female_fitness_contrast_X >0& male_fitness_contrast_Y >0| contrast =="X (female) vs Z (female)"& female_fitness_contrast_X >0& female_fitness_contrast_Z >0| contrast =="X (female) vs Z (male)"& female_fitness_contrast_X >0& male_fitness_contrast_Z >0| contrast =="X (female) vs W or Cytoplasmic (female)"& female_fitness_contrast_X >0&`female_fitness_contrast_W or Cytoplasmic`>0| contrast =="X (female) vs Cytoplasmic (male)"& female_fitness_contrast_X >0&`male_fitness_contrast_W or Cytoplasmic`>0| contrast =="Z (female) vs Autosome (male)"& female_fitness_contrast_Z >0& male_fitness_contrast_Autosome >0| contrast =="Z (female) vs X (male)"& female_fitness_contrast_Z >0& male_fitness_contrast_X >0| contrast =="Z (female) vs Y (male)"& female_fitness_contrast_Z >0& male_fitness_contrast_Y >0| contrast =="Z (female) vs Z (male)"& female_fitness_contrast_Z >0& male_fitness_contrast_Z >0| contrast =="Z (female) vs W or Cytoplasmic (female)"& female_fitness_contrast_Z >0&`female_fitness_contrast_W or Cytoplasmic`>0| contrast =="Z (female) vs Cytoplasmic (male)"& female_fitness_contrast_Z >0&`male_fitness_contrast_W or Cytoplasmic`>0| contrast =="W or Cytoplasmic (female) vs Autosome (male)"&`female_fitness_contrast_W or Cytoplasmic`>0& male_fitness_contrast_Autosome >0| contrast =="W or Cytoplasmic (female) vs X (male)"&`female_fitness_contrast_W or Cytoplasmic`>0& male_fitness_contrast_X >0| contrast =="W or Cytoplasmic (female) vs Y (male)"&`female_fitness_contrast_W or Cytoplasmic`>0& male_fitness_contrast_Y >0| contrast =="W or Cytoplasmic (female) vs Z (male)"&`female_fitness_contrast_W or Cytoplasmic`>0& male_fitness_contrast_Z >0| contrast =="W or Cytoplasmic (female) vs Cytoplasmic (male)"&`female_fitness_contrast_W or Cytoplasmic`>0&`male_fitness_contrast_W or Cytoplasmic`>0| contrast =="Autosome (male) vs X (male)"& male_fitness_contrast_Autosome >0& male_fitness_contrast_X >0| contrast =="Autosome (male) vs Y (male)"& male_fitness_contrast_Autosome >0& male_fitness_contrast_Y >0| contrast =="Autosome (male) vs Z (male)"& male_fitness_contrast_Autosome >0& male_fitness_contrast_Z >0| contrast =="Autosome (male) vs Cytoplasmic (male)"& male_fitness_contrast_Autosome >0&`male_fitness_contrast_W or Cytoplasmic`>0| contrast =="X (male) vs Y (male)"& male_fitness_contrast_X >0& male_fitness_contrast_Y >0| contrast =="X (male) vs Z (male)"& male_fitness_contrast_X >0& male_fitness_contrast_Z >0| contrast =="X (male) vs Cytoplasmic (male)"& male_fitness_contrast_X >0&`male_fitness_contrast_W or Cytoplasmic`>0| contrast =="Y (male) vs Z (male)"& male_fitness_contrast_Y >0& male_fitness_contrast_Z >0| contrast =="Y (male) vs Cytoplasmic (male)"& male_fitness_contrast_Y >0&`male_fitness_contrast_W or Cytoplasmic`>0| contrast =="Z (male) vs Cytoplasmic (male)"& male_fitness_contrast_Z >0&`male_fitness_contrast_W or Cytoplasmic`>0~"Inbreeding favoured in both contexts",# deleterious in both contexts contrast =="Autosome (female) vs Autosome (male)"& female_fitness_contrast_Autosome <0& male_fitness_contrast_Autosome <0| contrast =="Autosome (female) vs X (female)"& female_fitness_contrast_Autosome <0& female_fitness_contrast_X <0| contrast =="Autosome (female) vs X (male)"& female_fitness_contrast_Autosome <0& male_fitness_contrast_X <0| contrast =="Autosome (female) vs Y (male)"& female_fitness_contrast_Autosome <0& male_fitness_contrast_Y <0| contrast =="Autosome (female) vs Z (female)"& female_fitness_contrast_Autosome <0& female_fitness_contrast_Z <0| contrast =="Autosome (female) vs Z (male)"& female_fitness_contrast_Autosome <0& male_fitness_contrast_Z <0| contrast =="Autosome (female) vs W or Cytoplasmic (female)"& female_fitness_contrast_Autosome <0&`female_fitness_contrast_W or Cytoplasmic`<0| contrast =="Autosome (female) vs Cytoplasmic (male)"& female_fitness_contrast_Autosome <0&`male_fitness_contrast_W or Cytoplasmic`<0| contrast =="X (female) vs Autosome (male)"& female_fitness_contrast_X <0& male_fitness_contrast_Autosome <0| contrast =="X (female) vs X (male)"& female_fitness_contrast_X <0& male_fitness_contrast_X <0| contrast =="X (female) vs Y (male)"& female_fitness_contrast_X <0& male_fitness_contrast_Y <0| contrast =="X (female) vs Z (female)"& female_fitness_contrast_X <0& female_fitness_contrast_Z <0| contrast =="X (female) vs Z (male)"& female_fitness_contrast_X <0& male_fitness_contrast_Z <0| contrast =="X (female) vs W or Cytoplasmic (female)"& female_fitness_contrast_X <0&`female_fitness_contrast_W or Cytoplasmic`<0| contrast =="X (female) vs Cytoplasmic (male)"& female_fitness_contrast_X <0&`male_fitness_contrast_W or Cytoplasmic`<0| contrast =="Z (female) vs Autosome (male)"& female_fitness_contrast_Z <0& male_fitness_contrast_Autosome <0| contrast =="Z (female) vs X (male)"& female_fitness_contrast_Z <0& male_fitness_contrast_X <0| contrast =="Z (female) vs Y (male)"& female_fitness_contrast_Z <0& male_fitness_contrast_Y <0| contrast =="Z (female) vs Z (male)"& female_fitness_contrast_Z <0& male_fitness_contrast_Z <0| contrast =="Z (female) vs W or Cytoplasmic (female)"& female_fitness_contrast_Z <0&`female_fitness_contrast_W or Cytoplasmic`<0| contrast =="Z (female) vs Cytoplasmic (male)"& female_fitness_contrast_Z <0&`male_fitness_contrast_W or Cytoplasmic`<0| contrast =="W or Cytoplasmic (female) vs Autosome (male)"&`female_fitness_contrast_W or Cytoplasmic`<0& male_fitness_contrast_Autosome <0| contrast =="W or Cytoplasmic (female) vs X (male)"&`female_fitness_contrast_W or Cytoplasmic`<0& male_fitness_contrast_X <0| contrast =="W or Cytoplasmic (female) vs Y (male)"&`female_fitness_contrast_W or Cytoplasmic`<0& male_fitness_contrast_Y <0| contrast =="W or Cytoplasmic (female) vs Z (male)"&`female_fitness_contrast_W or Cytoplasmic`<0& male_fitness_contrast_Z <0| contrast =="W or Cytoplasmic (female) vs Cytoplasmic (male)"&`female_fitness_contrast_W or Cytoplasmic`<0&`male_fitness_contrast_W or Cytoplasmic`<0| contrast =="Autosome (male) vs X (male)"& male_fitness_contrast_Autosome <0& male_fitness_contrast_X <0| contrast =="Autosome (male) vs Y (male)"& male_fitness_contrast_Autosome <0& male_fitness_contrast_Y <0| contrast =="Autosome (male) vs Z (male)"& male_fitness_contrast_Autosome <0& male_fitness_contrast_Z <0| contrast =="Autosome (male) vs Cytoplasmic (male)"& male_fitness_contrast_Autosome <0&`male_fitness_contrast_W or Cytoplasmic`<0| contrast =="X (male) vs Y (male)"& male_fitness_contrast_X <0& male_fitness_contrast_Y <0| contrast =="X (male) vs Z (male)"& male_fitness_contrast_X <0& male_fitness_contrast_Z <0| contrast =="X (male) vs Cytoplasmic (male)"& male_fitness_contrast_X <0&`male_fitness_contrast_W or Cytoplasmic`<0| contrast =="Y (male) vs Z (male)"& male_fitness_contrast_Y <0& male_fitness_contrast_Z <0| contrast =="Y (male) vs Cytoplasmic (male)"& male_fitness_contrast_Y <0&`male_fitness_contrast_W or Cytoplasmic`<0| contrast =="Z (male) vs Cytoplasmic (male)"& male_fitness_contrast_Z <0&`male_fitness_contrast_W or Cytoplasmic`<0~"Inbreeding deleterious in both contexts",.default ="Intragenomic conflict")) %>%mutate(Evolutionary_conflict =if_else(is.na(Evolutionary_conflict), 0, Evolutionary_conflict),contrast =fct_relevel(contrast, "W or Cytoplasmic (female) vs Cytoplasmic (male)","X (female) vs Cytoplasmic (male)", "X (female) vs W or Cytoplasmic (female)","Autosome (female) vs Cytoplasmic (male)","Autosome (female) vs W or Cytoplasmic (female)","Z (female) vs Cytoplasmic (male)","Z (female) vs W or Cytoplasmic (female)","X (male) vs Cytoplasmic (male)","W or Cytoplasmic (female) vs X (male)","Autosome (male) vs Cytoplasmic (male)","W or Cytoplasmic (female) vs Autosome (male)","Z (male) vs Cytoplasmic (male)","W or Cytoplasmic (female) vs Z (male)","Y (male) vs Cytoplasmic (male)","W or Cytoplasmic (female) vs Y (male)","Autosome (female) vs X (female)","X (female) vs Z (female)","Autosome (female) vs Z (female)","X (female) vs X (male)","X (female) vs Autosome (male)","Autosome (female) vs X (male)","X (female) vs Z (male)","Autosome (male) vs X (male)","Z (female) vs X (male)","Autosome (female) vs Autosome (male)","Autosome (female) vs Z (male)","Z (female) vs Autosome (male)","Z (female) vs Z (male)","X (male) vs Z (male)","Autosome (male) vs Z (male)","X (female) vs Y (male)","Autosome (female) vs Y (male)","Z (female) vs Y (male)","X (male) vs Y (male)","Autosome (male) vs Y (male)","Y (male) vs Z (male)" ))make_genomic_conflict_plot_2 <-function(data, colour_pal){ data %>%ggplot(aes(x = r, y = D)) +geom_blank() +geom_tile(data = data %>%filter(relationship =="Intragenomic conflict"),aes(fill = Evolutionary_conflict*-1)) +#geom_tile(data = data,# aes(fill = Evolutionary_conflict)) + scale_fill_gradientn(colours = colour_pal, limits =c(0, 0.6), #na.value = "white",breaks =c(0, 0.6),labels =c("Weaker conflict", "Stronger conflict")) +labs(fill ="Evolutionary conflict") +new_scale_fill() +geom_tile(data = data %>%filter(relationship !="Intragenomic conflict"),aes(fill = relationship), alpha =0.75) +scale_fill_manual(values =c("#fbe6c5", "#d2fbd4"), labels =c("Inbreeding deleterious\nin both contexts", "Inbreeding favoured\nin both contexts")) +stat_contour(aes(z = Evolutionary_conflict*-1), colour ="black", binwidth =25,breaks =c(0, 0.2, 0.4, 0.6)) +facet_wrap(~contrast, nrow =6,scales ="free", strip.position =c("top"),labeller =label_glue('{`contrast`}')) +labs(x ="Individual-level relatedness coefficient",y ="Inbreeding depression",fill ="Evolutionary concordance") +scale_x_continuous(expand =c(0, 0), breaks =c(0, 0.25, 0.5, 0.75)) +scale_y_continuous(expand =c(0, 0)) +# labels = c(0, 25, 50, 75, 90)) +theme(panel.border =element_rect(fill =NA, colour ="black", size = .8),panel.grid.minor =element_blank(),strip.background =element_rect(colour ="black", fill ="Aliceblue", linewidth = .8),strip.text =element_text(size =6),axis.title.x =element_markdown(),axis.text =element_text(size =8),#plot.title = element_text(hjust = 0.5, size = 8),legend.position ="none") }make_genomic_conflict_plot_2(Intragenomic_all_combos, oranges)
Make Figure 2
We can also plot the cost to female direct reproductive output:
Figure 2. direct female reproductive output for inbreeding females at the threshold where chromosome-specific selection should no longer favour inbreeding. Reproductive output is expressed relative to females that exclusively outcross. Cytoplasmic chromosomes are included here to illustrate that the inclusive fitness of alleles on these chromosomes is absolutely aligned with individual-level direct female fitness.
How fast do inbreeding families become homozygous?
Run a quick set of simulations to work this out, where female-male pairs of all possible genotypes at the inbreeding locus begin new families. We simulate the autosomal, X and Z cases, as these are where the number of I alleles carried by individuals matters.
A sampling function that can handle vectors of length one
Code
# so we can sample from vectors with length 1 without this being interpreted as an integersample_vec <-function(x, ...) x[sample(length(x), ...)]
A function that builds the appropriate inheritance system.
Code
make_mating_table <-function(gene_location){ make_offspring <-function(X, Y, offspring_genotype, zygote_freq, sex){data.frame(Female_genotype = X,Male_genotype = Y, offspring_genotype, zygote_freq, sex) }# Specify the possible offspring genotypes for all the potential crosses; we use these for the offspring_genotype argument in the make_offspring function# offspring genotypes offspring_genotypes_1 <-c(2,2) offspring_genotypes_2 <-c(1, 1, 2, 2) offspring_genotypes_3 <-c(1, 1) offspring_genotypes_4 <-c(0, 0, 1, 1, 2, 2) offspring_genotypes_5 <-c(0, 0, 1, 1) offspring_genotypes_6 <-c(0, 0) offspring_genotypes_7 <-c(1, 2) offspring_genotypes_8 <-c(0, 1, 1, 2) offspring_genotypes_9 <-c(0, 1) offspring_genotypes_10 <-c(1, 0) # this is diff from above bc of the order with the sexes offspring_genotypes_11 <-c(2, 1) offspring_genotypes_12 <-c(2,1,1,0)# offspring sex offspring_sex_2 <-c(0, 1) offspring_sex_4 <-c(0, 1, 0, 1) offspring_sex_6 <-c(0, 1, 0, 1, 0, 1)# even frequency of two offspring genotypes freq_2 <-rep(0.5, 2)# even frequency between four offspring types freq_4 <-rep(0.25, 4)# when there are 6 offspring genotypes freq_6 <-c(0.125, 0.125,0.25, 0.25,0.125, 0.125)if(gene_location =="A"){ books <-rbind(make_offspring(2, 2, offspring_genotypes_1, freq_2, offspring_sex_2),make_offspring(2, 1, offspring_genotypes_2, freq_4, offspring_sex_4),make_offspring(2, 0, offspring_genotypes_3, freq_2, offspring_sex_2),make_offspring(1, 2, offspring_genotypes_2, freq_4, offspring_sex_4),make_offspring(1, 1, offspring_genotypes_4, freq_6, offspring_sex_6),make_offspring(1, 0, offspring_genotypes_5, freq_4, offspring_sex_4),make_offspring(0, 2, offspring_genotypes_3, freq_2, offspring_sex_2),make_offspring(0, 1, offspring_genotypes_5, freq_4, offspring_sex_4),make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2) ) }if(gene_location =="X"){ books <-rbind(make_offspring(2, 1, offspring_genotypes_7, freq_2, offspring_sex_2),make_offspring(2, 0, offspring_genotypes_3, freq_2, offspring_sex_2),make_offspring(1, 1, offspring_genotypes_8, freq_4, offspring_sex_4),make_offspring(1, 0, offspring_genotypes_5, freq_4, offspring_sex_4),make_offspring(0, 1, offspring_genotypes_9, freq_2, offspring_sex_2),make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2) ) }if(gene_location =="Y"){ books <-rbind(make_offspring(0, 1, offspring_genotypes_10, freq_2, offspring_sex_2),make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2) ) }if(gene_location =="Z"){ books <-rbind(make_offspring(1, 2, offspring_genotypes_11, freq_2, offspring_sex_2),make_offspring(0, 2, offspring_genotypes_3, freq_2, offspring_sex_2),make_offspring(1, 1, offspring_genotypes_12, freq_4, offspring_sex_4),make_offspring(0, 1, offspring_genotypes_5, freq_4, offspring_sex_4),make_offspring(1, 0, offspring_genotypes_10, freq_2, offspring_sex_2),make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2) ) }if(gene_location =="W"){ books <-rbind(make_offspring(1, 0, offspring_genotypes_9, freq_2, offspring_sex_2),make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2) ) }if(gene_location =="C"){ books <-rbind(make_offspring(1, 0, offspring_genotypes_3, freq_2, offspring_sex_2),make_offspring(1, 1, offspring_genotypes_3, freq_2, offspring_sex_2),make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2),make_offspring(0, 1, offspring_genotypes_6, freq_2, offspring_sex_2) ) }if(gene_location =="P"){ books <-rbind(make_offspring(1, 0, offspring_genotypes_6, freq_2, offspring_sex_2),make_offspring(1, 1, offspring_genotypes_3, freq_2, offspring_sex_2),make_offspring(0, 0, offspring_genotypes_6, freq_2, offspring_sex_2),make_offspring(0, 1, offspring_genotypes_3, freq_2, offspring_sex_2) ) }return(books) }offspring_genotypes_autosome <-make_mating_table("A")offspring_genotypes_X <-make_mating_table("X")offspring_genotypes_Y <-make_mating_table("Y")offspring_genotypes_Z <-make_mating_table("Z")offspring_genotypes_W <-make_mating_table("W")offspring_genotypes_C <-make_mating_table("C")offspring_genotypes_P <-make_mating_table("P")
A function that takes two parental genotypes and produces offspring
Code
sample_mating_table <-function(inheritance_scheme, f, mother){# cut to possible genotypes possibilities <- inheritance_scheme[inheritance_scheme$Female_genotype == mother[4] & inheritance_scheme$Male_genotype == mother[9], c(3,5)]# get prob of producing each genotype probs <- inheritance_scheme[inheritance_scheme$Female_genotype == mother[4] & inheritance_scheme$Male_genotype == mother[9], 4]# sample possibilities[sample(size = f,x =nrow(possibilities), prob = probs,replace =TRUE), ]}
Load the parameter space
Code
resolution <-25starting_pop_size_autosomes <-2000# both sexes harbour two copies of each autosomal chromosome = 1000 autosomal haplotypesparameters <-expand_grid(chromosome =c("A", "X", "Y", "Z", "W", "C", "P"),v =c(8, 80),D =seq(0, -0.99, length = resolution), # inbreeding depressionrefractory_period_prop_cohort_alive =seq(0.01, 1, length = resolution) ) %>%full_join(tibble(chromosome =c("A", "A", "A", "A", "A", "A","X", "X", "X", "X","Y","Z", "Z", "Z", "Z","W","C", "C","P", "P"),sex_expressed =c(0, 0, 0, 1, 1, 1,0, 1, 1, 1,0,0, 0, 0, 1,1,0, 1,0, 1),dominance =c(0, 0.5, 1, 0, 0.5, 1,1, 0, 0.5, 1, 1, 0, 0.5, 1, 1,1,1, 1,1, 1)) %>%mutate(Starting_pop_size =case_when(chromosome =="A"~ starting_pop_size_autosomes, chromosome =="X"| chromosome =="Z"~ starting_pop_size_autosomes /0.75, chromosome =="Y"| chromosome =="W"~ starting_pop_size_autosomes*4, chromosome =="C"| chromosome =="P"~ starting_pop_size_autosomes*2)),relationship ="many-to-many", by ="chromosome") %>%mutate(baseline_mean_lifespan =1,v = v / (Starting_pop_size /2),f =5, refractory_period =-log(refractory_period_prop_cohort_alive),mutation_time =5, # this is when the mutation can be introduced fromtime_end =1000, # with avg lifespan = 1, this is ~ roughly 1000 gensparameter_space_ID =row_number(),mutation_events =5)parameters_autosome <- parameters %>%filter(chromosome =="A") %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_X <- parameters %>%filter(chromosome =="X") %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_Y <- parameters %>%filter(chromosome =="Y") %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_Z <- parameters %>%filter(chromosome =="Z") %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_W <- parameters %>%filter(chromosome =="W") %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_C <- parameters %>%filter(chromosome =="C") %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_P <- parameters %>%filter(chromosome =="P") %>%slice_sample(prop =1) # shuffle to equalise workload across jobsif(!file.exists("parameters/parameters_autosome.txt")){ parameters_autosome %>%write.table("parameters/parameters_autosome.txt") parameters_X %>%write.table("parameters/parameters_X.txt") parameters_Y %>%write.table("parameters/parameters_Y.txt") parameters_Z %>%write.table("parameters/parameters_Z.txt") parameters_W %>%write.table("parameters/parameters_W.txt") parameters_C %>%write.table("parameters/parameters_C.txt") parameters_P %>%write.table("parameters/parameters_P.txt")}
The simulation function
Code
continuous_time_simulation <-function(row, parameters, inheritance_scheme){#print(paste("Doing row", row)) # this shows which row in the parameter space is being modelled Starting_pop_size <-round(parameters$Starting_pop_size[row], 0) f <- parameters$f[row] # fecundity constant mutation_time <- parameters$mutation_time[row] # introduce an I allele after family structure is established baseline_mean_lifespan <- parameters$baseline_mean_lifespan[row] # constant at 1 time_end <- parameters$time_end[row] # a cut-off point for each run sex_expressed <- parameters$sex_expressed[row] chromosome <- parameters$chromosome[row] v <- parameters$v[row] refractory_period <- parameters$refractory_period[row] D <- parameters$D[row] dominance <- parameters$dominance[row] parameter_space_ID <- parameters$parameter_space_ID[row] mutation_events <- parameters$mutation_events[row]# Set the number of breeding sites breeding_sites <-round(0.2*Starting_pop_size, 0)# what inheritance system does this run follow offspring_genotypes <- inheritance_scheme# Set the maximum number of I alleles that can be found in each sexif(chromosome =="A"){ female_max_I <-2 male_max_I <-2 }if(chromosome =="X"){ female_max_I <-2 male_max_I <-1 }if(chromosome =="Z"){ female_max_I <-1 male_max_I <-2 }if(chromosome =="Y"){ female_max_I <-0 male_max_I <-1 }if(chromosome =="W"){ female_max_I <-1 male_max_I <-0 }if(chromosome =="C"| chromosome =="P"){ female_max_I <-1 male_max_I <-1 }# make matrix to hold results; updated as sim progresses# col1 = time, col2 = prop I, col3 = pop size, col4 = prop virgin female deaths results_matrix <-matrix(nrow = time_end*4+2, ncol =4) # record each time point# make matrix to hold population; updated as sim progresses# col1 = ID # col2 = Family ID# col3 = Sex: females = 1 and males = 0# col4 = Genotype: 0, 1 and 2 = copies of inbreeding allele# col5 = mortality rate# col6 = encountered relative: NA = NO, 1 = YES# col7 = mating state: -Inf not in pop, NA = unmated, real = out, Inf = mated female# col8 = inbred mating: NA = NO, 1 = YES (only matters for females)# col9 = mated_genotype: NA = unmated, otherwise 0,1,2 (see mating table)# col10 = breeding site: NA = NO, 1 = YES # col11 = no. matings (only matters for males)# col12 = offspring produced: NA = NO, 1 = YES pop_matrix <-matrix(nrow = Starting_pop_size*2, # pop expands with initial repro pulsencol =12)# ID & Family ID pop_matrix[1:Starting_pop_size, 1:2] <-1:Starting_pop_size# assign sex pop_matrix[1:Starting_pop_size, 3] <-rbinom(n = Starting_pop_size, 1, prob =0.5)# female_starting_genotype pop_matrix[pop_matrix[,3] <1&!is.na(pop_matrix[,3]), 4] <-0# male_starting_genotype pop_matrix[pop_matrix[,3] >0&!is.na(pop_matrix[,3]), 4] <-0# assign mortality rates pop_matrix[1:Starting_pop_size, 5] <-1/baseline_mean_lifespan# set the unused rows to state -Inf pop_matrix[(Starting_pop_size +1):nrow(pop_matrix), 7] <--Inf# mate count pop_matrix[1:Starting_pop_size, 11] <-0# offspring production status pop_matrix[1:Starting_pop_size, 12] <-NA# populate breeding sites# the starting no. of females generally exceeds the number of breeding sites, which is starting_pop_size/f. The code below selects the initial breeding site holdersif(nrow(pop_matrix[pop_matrix[,3] >0&!is.na(pop_matrix[,3]),]) > breeding_sites){ initial_breeders <-head(pop_matrix[pop_matrix[,3] >0&!is.na(pop_matrix[,3]),1], breeding_sites) pop_matrix[initial_breeders,10] <-1# take advantage of ID = row number for initial pop } else{pop_matrix[pop_matrix[,3] >0&!is.na(pop_matrix[,3]), 10] <-1}# Initialise counter for the results table next_update <-0# keep track of when to update the results next_row <-0# keep track of which row to update# Initialise the Individual_ID and Family_ID counters Individual_ID_counter <- Starting_pop_size Family_ID_counter <- Starting_pop_size # each individual descends from a distinct family at onset# Initialise the timer t t <-0# Set initial pop size and freq of I allele for results table Prop_I <-0 pop_size <- Starting_pop_size total_female_deaths <-0 mated_female_deaths <-0# Start population without the I allele to generate family structure# Flips to 1 at mutant intro time point mutant_introduced <-0 keep_going <-TRUE# if the inbreeding allele fixes or goes extinct, this will change to false and the while loop will quit early# With the initial population ready to go, start the timer and let the simulation run.while(t <= time_end & keep_going){print(paste0("Population size = ", pop_size, ", breeders = ", sum(pop_matrix[,10], na.rm = T), ", time = ", round(t, 3), ", Prop I =", Prop_I, ", mutation events =", mutant_introduced))# find next event # next death: this is the sum of the mortality rates for all individuals in the population next_death <- t +rexp(n =1, rate =sum(pop_matrix[, 5], na.rm = T))# next receptive mating encounter# find no. of females in mating pool & separate by encounter experience receptive_females_first_encounter <- pop_matrix[pop_matrix[,3] >0&is.na(pop_matrix[,6]) &is.na(pop_matrix[,7]),, drop =FALSE] receptive_females_second_encounter <- pop_matrix[pop_matrix[,3] >0&!is.na(pop_matrix[,6]) &is.na(pop_matrix[,7]),, drop =FALSE]# find no. of males in mating pool receptive_males <- pop_matrix[pop_matrix[,3] <1&is.na(pop_matrix[,7]),, drop =FALSE]# Find the time the next encounter occurs: plug the sum of the rates into the exponential function. # The population level encounter rate is the product of the rate at which a single male finds a single female, the number of receptive females in the population, and the number of receptive males in the population next_first_encounter <- t +rexp(n =1, rate = v*nrow(receptive_females_first_encounter)*nrow(receptive_males)) next_secondary_encounter <- t +rexp(n =1, rate = v*nrow(receptive_females_second_encounter)*nrow(receptive_males))# time in - Inf, Inf and Na are possible options that the code can handle next_time_in <-min(pop_matrix[is.finite(pop_matrix[,7]),7])# find which event happens next and update t t <-pmin(next_death, next_time_in, next_first_encounter, next_secondary_encounter, next_update, # update the populationna.rm =TRUE) # ... if a rate is 0, NaN produced.if(t == next_update &!is.na(next_update)){# record time, I prop and pop size results_matrix[next_row+1,1] <- t results_matrix[next_row+1,2] <-round(Prop_I, 4) results_matrix[next_row+1,3] <- pop_size # popsize results_matrix[next_row+1,4] <-round(mated_female_deaths / total_female_deaths, 3) next_update <- next_update +0.25 next_row <- next_row +1 total_female_deaths <-0# reset the count mated_female_deaths <-0# reset the count }if(t == next_death){# remove an individual from the pop who_died <-sample_vec(size =1, # choose onex = pop_matrix[!is.na(pop_matrix[,1]),1], # subset to current popprob = pop_matrix[!is.na(pop_matrix[,5]),5]) # weight by mortality rate# add a death if it was a femaleif(nrow(pop_matrix[pop_matrix[,1] == who_died &!is.na(pop_matrix[,1]) & pop_matrix[,3] >0,, drop =FALSE]) >0){total_female_deaths <- total_female_deaths +1}# add virgin female deathsif(nrow(pop_matrix[pop_matrix[,1] == who_died &!is.na(pop_matrix[,1]) & pop_matrix[,3] >0&is.infinite(pop_matrix[,7]),, drop =FALSE]) >0){mated_female_deaths <- mated_female_deaths +1}# remove individual from pop matrix pop_matrix[pop_matrix[,1] == who_died, 7] <--Inf# NA means time-in here, so special edit required pop_matrix[pop_matrix[,1] == who_died, c(1:6, 8:12)] <-NA# re-order to make steps like adding offspring easier later on pop_matrix <- pop_matrix[order(pop_matrix[,1]),] }# check if there are free breeding sites and whether females are available to fill them current_breeders <-sum(pop_matrix[, 10], na.rm = T)# get list of IDs for floating females floating_females <- pop_matrix[!is.na(pop_matrix[,1]) &# alive pop_matrix[,3] >0&# femaleis.na(pop_matrix[,10]), # non-breeding1] # return the IDs only# If so, recruit a new breeder# All prospective females have equal probabilityif(current_breeders < breeding_sites &length(floating_females) >0){# assign the new breeder new_breeder <-sample_vec(size =1, # choose onex = floating_females) # subset to floaters pop_matrix[pop_matrix[,1] == new_breeder, 10] <-1 }if(t == next_time_in &!is.na(next_time_in)){ # a male re-enters the mating pool pop_matrix[pop_matrix[,7] == next_time_in, 7] <-NA# change to receptive }#### matingif(t == next_first_encounter &!is.na(next_first_encounter)){# does first encounter lead to (inbred) mating?# Determine whether a heterozygote inbreeds on this occasion. # Depends on genotype if this matters heterozygote_inbreeds <-rbinom(1, 1, prob = dominance)# which female female_ID <-sample_vec(receptive_females_first_encounter[,1], 1)# get meta-data female <-subset(pop_matrix, pop_matrix[,1] == female_ID)# how many inbreeding alleles does she carry? alleles_female <- female[,4] mates <-NULL# reset this every time as a safeguard - MAYBE REMOVE?# find brothers that are in the mating pool brothers <- pop_matrix[pop_matrix[,2] == female[, 2] &# find family members pop_matrix[,3] <1&# that are maleis.na(pop_matrix[,7]) &# and in the mating pool!is.na(pop_matrix[,1]), # remove NAs1] # find the specific brother - if there aren't any, inbreeding does not happenif(length(brothers) >0){# choose brother randomly chosen_brother <-subset(pop_matrix, pop_matrix[,1] ==sample_vec(size =1, x = brothers))# how many inbreeding alleles does he carry? alleles_brother <- chosen_brother[,4] brother_ID <- chosen_brother[,1] }else{alleles_brother <-0} # we need this for the next if statement# now determine whether inbreeding occurs:# which individual expresses the allele# does that individual have the allele# is it expressed (depends on genomic region, no. copies and dominance)if(# female expression determines outcome# dominance doesn't matterlength(brothers) >0& sex_expressed >0& female_max_I == alleles_female |# dominance matterslength(brothers) >0& sex_expressed >0&0< alleles_female & alleles_female < female_max_I & heterozygote_inbreeds >0|# male expression determines outcome# dominance doesn't matterlength(brothers) >0& sex_expressed <1& male_max_I == alleles_brother |# dominance matterslength(brothers) >0& sex_expressed <1&0< alleles_brother & alleles_brother < male_max_I & heterozygote_inbreeds >0){# do inbreeding# update the pop matrix# female pop_matrix[pop_matrix[,1] == female_ID, 6] <-1# relative has been encountered pop_matrix[pop_matrix[,1] == female_ID, 7] <-Inf# female leaves mating pool pop_matrix[pop_matrix[,1] == female_ID, 8] <-1# inbreeding occurs pop_matrix[pop_matrix[,1] == female_ID, 9] <- alleles_brother # mates genotype# male pop_matrix[pop_matrix[,1] == brother_ID, 7] <- t + refractory_period # male leaves mating pool pop_matrix[pop_matrix[,1] == brother_ID, 8] <-1# inbreeding occurs pop_matrix[pop_matrix[,1] == brother_ID &!is.na(pop_matrix[,1]), 11] <- pop_matrix[pop_matrix[,1] == brother_ID &!is.na(pop_matrix[,1]), 11] +1 } else{# inbreeding is avoided# females that had no receptive brother to encounter are recorded as having had their chance for inbreeding early in life. When the male refractory period != 0, this is possible but unlikely (because all siblings are produced at the same time). Most commonly, this will occur when a female produces an all-female brood (0.03125 probability when f=5) pop_matrix[pop_matrix[,1] == female_ID, 6] <-1# relative has been encountered } }if(t == next_secondary_encounter &!is.na(next_secondary_encounter)){ # If the individual has already encountered a sibling, don't swap and let encounter proceed. # which female female_ID <-sample_vec(receptive_females_second_encounter[,1], 1)# get meta-data female <-subset(pop_matrix, pop_matrix[,1] == female_ID)# how many inbreeding alleles does she carry? alleles_female <- female[,4] # which male male_ID <-sample_vec(receptive_males[,1], 1)# get meta-data male <-subset(pop_matrix, pop_matrix[,1] == male_ID)# how many inbreeding alleles does he carry? alleles_male <- male[,4] # If the pair happen to be siblings, check if they inbreed # Determine whether a heterozygote inbreeds on this occasion. # Depends on genotype if this matters heterozygote_inbreeds <-rbinom(1, 1, prob = dominance)if(# female expression determines outcome# dominance doesn't matter female[,2] == male[,2] & sex_expressed >0& female_max_I == alleles_female |# dominance matters female[,2] == male[,2] & sex_expressed >0&0< alleles_female & alleles_female < female_max_I & heterozygote_inbreeds >0|# male expression determines outcome# dominance doesn't matter female[,2] == male[,2] & sex_expressed <1& male_max_I == alleles_male |# dominance matters female[,2] == male[,2] & sex_expressed <1&0< alleles_male & alleles_male < male_max_I & heterozygote_inbreeds >0){# do inbreeding# update the pop matrix# female pop_matrix[pop_matrix[,1] == female_ID, 7] <-Inf# female leaves mating pool pop_matrix[pop_matrix[,1] == female_ID, 8] <-1# inbreeding occurs pop_matrix[pop_matrix[,1] == female_ID, 9] <- alleles_male # mates genotype# male pop_matrix[pop_matrix[,1] == male_ID, 7] <- t + refractory_period # male leaves mating pool pop_matrix[pop_matrix[,1] == male_ID &!is.na(pop_matrix[,1]), 11] <- pop_matrix[pop_matrix[,1] == male_ID &!is.na(pop_matrix[,1]), 11] +1 } else{# do outbreeding# update the pop matrix# female pop_matrix[pop_matrix[,1] == female_ID, 7] <-Inf# female leaves mating pool pop_matrix[pop_matrix[,1] == female_ID, 9] <- alleles_male # mates genotype# male pop_matrix[pop_matrix[,1] == male_ID, 7] <- t + refractory_period # male leaves mating pool pop_matrix[pop_matrix[,1] == male_ID &!is.na(pop_matrix[,1]), 11] <- pop_matrix[pop_matrix[,1] == male_ID &!is.na(pop_matrix[,1]), 11] +1 } }# Consequences of death and mating: reproduction# check if a female can now produce offspring, either because they're previously mated and have secured a breeding site or because they already hold a breeding site and have now mated# make sure that previous breeders are excluded new_mated_breeder <- pop_matrix[is.infinite(pop_matrix[,7]) &# mated!is.na(pop_matrix[,10]) &# holds breeding siteis.na(pop_matrix[,12]),, drop =FALSE] # hasn't reproducedif(nrow(new_mated_breeder) >0){# add offspring to the population# each mated female that holds a breeding site produces f offspring# first check whether the mutant I allele should be addedif(mutant_introduced < mutation_events & t > mutation_time){ which_sex <-rbinom(1, 1, prob =0.5)if(chromosome =="A"& which_sex ==1| chromosome =="X"& which_sex ==1| chromosome =="Z"& which_sex ==1){ new_mated_breeder[4] <-1 }if(chromosome =="A"& which_sex ==0| chromosome =="X"& which_sex ==0| chromosome =="Z"& which_sex ==0){ new_mated_breeder[9] <-1 }if(chromosome =="W"| chromosome =="C"){ new_mated_breeder[4] <-1 }if(chromosome =="Y"| chromosome =="P"){ new_mated_breeder[9] <-1 } mutant_introduced <- mutant_introduced +1 } next_row_to_fill <-length(pop_matrix[!is.na(pop_matrix[,1]),1]) +1 last_row_to_fill <- next_row_to_fill + f -1 next_ID <- Individual_ID_counter +1 last_ID <- Individual_ID_counter + f Family_ID_counter <- Family_ID_counter +1# assign IDs pop_matrix[next_row_to_fill:last_row_to_fill, 1] <- next_ID:last_ID# assign all offspring to a single family pop_matrix[next_row_to_fill:last_row_to_fill, 2] <- Family_ID_counter# assign sex and genotype using our mating table sampling function offspring_genos <-sample_mating_table(inheritance_scheme, f, mother = new_mated_breeder) pop_matrix[next_row_to_fill:last_row_to_fill, 3] <- offspring_genos[,2] pop_matrix[next_row_to_fill:last_row_to_fill, 4] <- offspring_genos[,1]# assign mortality ratesif(is.na(new_mated_breeder[8])){ pop_matrix[next_row_to_fill:last_row_to_fill, 5] <-1/baseline_mean_lifespan } else{ # apply effect of inbreeding depression pop_matrix[next_row_to_fill:last_row_to_fill, 5] <-1/(baseline_mean_lifespan + D) }# fill in the mating and breeding site details - everyone starts as a floating virgin pop_matrix[next_row_to_fill:last_row_to_fill, 6:10] <-NA# mate count pop_matrix[next_row_to_fill:last_row_to_fill, 11] <-0# update the mothers offspring production status pop_matrix[pop_matrix[,1] == new_mated_breeder[1], 12] <-1# update the individual ID counter (redundant but more readable to do this here) Individual_ID_counter <- last_ID } # Calculate the frequency of the I allele, quit early if I fixes or goes extinct pop_size <-nrow(pop_matrix[!is.na(pop_matrix[,1]),, drop =FALSE]) # use this to update the results n_females <-nrow(pop_matrix[!is.na(pop_matrix[,1]) & pop_matrix[,3] >0,, drop =FALSE]) n_males <- pop_size - n_females# calc allele freq if autosomal locus if(chromosome =="A"){ Prop_I <-sum(pop_matrix[,4], na.rm = T)/(pop_size*2) # x2 because diploid }# calc allele freq if W locus if(chromosome =="W"){ Prop_I <-sum(pop_matrix[,4], na.rm = T)/n_females }# calc allele freq if Y locus if(chromosome =="Y"){ Prop_I <-sum(pop_matrix[,4], na.rm = T)/n_males }# calc allele freq if X locus if(chromosome =="X"){ Prop_I <-sum(pop_matrix[,4], na.rm = T)/(n_females*2+ n_males) }# calc allele freq if Z locus if(chromosome =="Z"){ Prop_I <-sum(pop_matrix[,4], na.rm = T)/(n_females + n_males*2) }# calc allele freq if C locus if(chromosome =="C"| chromosome =="P"){ Prop_I <-sum(pop_matrix[,4], na.rm = T)/pop_size }# quit conditionif(mutant_introduced >0& Prop_I >0.9| mutant_introduced >0& Prop_I ==0| pop_size <2){keep_going <-FALSE} } results_matrix[next_row+1,1] <- t results_matrix[next_row+1,2] <-round(Prop_I, 4) results_matrix[next_row+1,3] <- pop_size results_matrix[next_row+1,4] <-round(mated_female_deaths / total_female_deaths, 3) results_matrix <- results_matrix[-(next_row+2:nrow(results_matrix)),]# save results as a csv. results_matrixwrite.csv(results_matrix,paste("results/rowID_", parameter_space_ID, chromosome, ".csv", sep =""))#write.csv(results_matrix,# paste("sim_results/rowID_", # parameter_space_ID, # chromosome, ".csv", # sep = ""))}
Run the simulation
In practice, I ran the simulations on JGU’s Mogon computing cluster. See the HPC_inbreeding_script.R and the batch script run_inbreeding_sim. To run the simulation for a single parameter space, you could run continuous_time_simulation(1, parameters_P, offspring_genotypes_P)
